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^ |Abstract 

5 Functional connectomes capture brain interactions via synchronized fluctuations in the functional magnetic resonance 

^ imaging signal. If measured during rest, they map the intrinsic functional architecture of the brain. With task-driven 

h* experiments they represent integration mechanisms between specialized brain areas. Analyzing their variability across 

"'subjects and conditions can reveal markers of brain pathologies and mechanisms underlying cognition. Methods of 

estimating functional connectomes from the imaging signal have undergone rapid developments and the literature is full 

" of diverse strategies for comparing them. This review aims to clarify links across functional-connectivity methods as well 

^ ! as to expose different steps to perform a group study of functional connectomes. 



< 



d 

•i-H 

en 



Keywords: Functional connectivity, connectome, group study, effective connectivity, fMRI, resting-state 



> 

o 

oo 

00 
CO 

o 

CO 



t 



1. Introduction 

Functional connectivity reveals the synchronization of 
distant neural systems via correlations in neurophysiolog- 
ical measures of brain activity [14, 37]. Given that high- 
level function emerges from the interaction of specialized 
units [110], functional connectivity is an essential part of 
the description of brain function, that complements the 
localizationist picture emerging from the systematic map- 
ping of regions recruited in tasks [101]. However, while 
there exists a well-defined standard analysis framework 
for activation mapping that enables statistically-controlled 
comparisons across subjects [39], group- level analysis of 
functional connectivity still face many open methodolog- 
ical challenges. Deriving a picture of a single subject's 
functional connectivity is by itself not straightforward, as 
the brain comprises a myriad of interacting subsystems 
and its connectivity must be decomposed into simplified 
and synthetic representations. An important view of brain 
connectivity is that of distributed functional networks de- 
picted by their spatial maps [31]. Another no less impor- 
tant and complementary view is that of connections link- 
ing localized functional modules depicted as a graph [17]. 
This representation of brain connectivity is often called 
the functional connectome [102] and is the focus of intense 
worldwide research efforts as it holds promises of new in- 
sights in cognition and pathologies [13, 30, 45]. 

The purpose of this paper is to review methodological 
progress in the estimation of functional connectomes from 
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blood oxygenation level dependent (BOLD) based func- 
tional magnetic resonance imaging (fMRI) data and their 
comparisons across individuals. It does not attempt to be 
exhaustive, as the field is wide and moving rapidly, but de- 
tails specific tools and guidelines that, in the experience of 
the authors, lead to controlled and powerful inter-subject 
comparisons. The paper is focused on functional connec- 
tomes in contrast to structural connectomes, as the infer- 
ence of functional connectivity requires important statisti- 
cal modeling considerations that are vastly different from 
the complications involved with estimating structural con- 
nectivity. While the notion of functional connectomics is 
often associated with the study of resting state [13], the 
methods presented in this paper are also relevant for task- 
based studies. On the other hand, the paper has a focus on 
fMRI; although the core concepts presented can be applied 
to magnetoencephalography (MEG) or electroencephalog- 
raphy (EEG) [103], additional specific problems such as 
source reconstruction must be considered [93]. 

"Functional connectivity" is defined as a measure of 
synchronization in brain signals [35]. More generally, it 
is interesting as a window on underlying synchrony on 
neural processes [63]. By "functional connectome", here 
we specifically denote a graph representing functional in- 
teractions in the brain, where the term "graph" is taken 
in its mathematical sense: a set of nodes connected to- 
gether by edges. Graph nodes (brain regions) correspond 
to spatially-contiguous and functionally-coherent patches 
of gray matter and edges describe long-range synchroniza- 
tions between nodes that are putatively subtended by large 
fiber pathways [68]. A graph can be weighted or not, 
and is completely equivalent to its adjacency matrix, a 
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symmetric matrix tabulating the connection weights be- 
tween each pair of nodes. Functional-connectivity graphs 
are used to represent evoked activity, as in task-response 
studies [72] , as well as ongoing activity, present in the ab- 
sence of specific tasks or in the background during task 
and often studied in so-called resting state experiments 
[83] . Another important notion that arises from the study 
of distributed modes of brain function is that of specialized 
functional networks 1 [31]. With our definition of the func- 
tional connectome, functional networks are not directly 
building blocks of the connectome but appear as a conse- 
quence of the graphical structure [116, 117]. 

The paper is organized as follows. First we discuss es- 
timation of functional connectomes. This part, akin to a 
first-level analysis in standard activation mapping method- 
ology, is not in itself a group-level operation, but it is a 
critical step for inter-subject comparison. In a following 
section, we discuss several strategies for comparing con- 
nectomes across subjects. Finally we discuss the links be- 
tween the representation of brain connectivity as graphs 
of functional connectivity and more complex models, such 
as effective-connectivity models. 

2. Estimating functional connectomes 

Here we discuss the inference of connectomes from 
functional brain imaging data. We start with preprocess- 
ing considerations, followed by the choice of nodes i.e. re- 
gions, signal extraction, and the estimation of graphs. 

2.1. Preprocessing considerations 

In addition to standard preprocessing performed 
for task-based analysis (slice-timing correction, realign- 
ment, spatial normalization, and possibly smoothing), 
connectivity-based analysis require additional denoising to 
separate intrinsic activity from confounding signals. This 
process involves regressing time series capturing sources 
of structured noise from the fMRI data. Physiological 
noise due to cardiac and respiration are two important 
noise signals [11, 12, 53, 67] that are difficult to control 
for and as a result are not commonly regressed out. In- 
stead the mean signal from white matter (WM) and cere- 
brospinal fluid (CSF) are used as surrogates to measure 
these sources of noise as well as other scanner induced sig- 
nal fluctuations [31, 67]. More complex models account for 
spatial variation in noise by incorporating voxel-specific re- 
gressors of neighboring WM (ANATICOR [55]) or the top 
components from a principal components analysis of high- 
variance signals (CompCor [7]). Head motion induced sig- 
nal fluctuations are accounted for by incorporating move- 
ment parameters [31, 41, 67]. The global mean time series 



has been proposed as an additional noise regressor that 
appears to improve the spatial specificity of connectivity 
results [31, 32]. This practice has become controversial 
since the global signal regression introduces negative corre- 
lations [19, 77, 90]. Removing these sources of nuisance in 
addition to linear trends results in more contrasted corre- 
lation matrices that improve the delineation of functional 
structures (fig. 2). 

Filtering to remove high frequencies is often performed, 
based on the initial observation that fluctuations impli- 
cated in resting-state functional connectivity are predom- 
inately slower than 0.1 Hz [14, 23]. While high-pass and 
low-pass filtering decrease the impact of some confounds, 
recent studies have shown that connectivity is present 
across the full spectrum of observed frequencies [99, 113]. 
Regressing out a good choice of confound signals is more 
specific than frequency filtering, and in our experience 
gives more contrasted correlation matrices 2 . In addition, 
the recent developments of very rapid acquisition proto- 
cols prevent aliasing of the physiological noise with the 
neural signal and give access to more specific noise con- 
founds than traditional low-TR sequences [16]. 

It is important to keep in mind that the proposed cor- 
rection strategies are approximate and not definitive tech- 
niques. This has become particularly apparent for head 
motion with reports that micromovements on the scale of 
< 0.2 mm can induce artefactual group- level findings even 
when motion is accounted for in preprocessing [82, 92, 112]. 
Special care must be taken to adequately control for resid- 
ual impact of head motion in the group model [92, 112]. 

2.2. Defining regions 

The choice of regions of interests (ROIs) that define 
the nodes of the graphs can be very important both in 
the estimation of connectomes and for group comparison 
[119]. Unsurprisingly, simulations have shown that ex- 
tracting signal from ROIs that did not match functional 
units would lead to erronous graph estimation [100]. Dif- 
ferent strategies to define suitable ROIs coexist. While 
dense parcellation approaches cover a large fraction of the 
brain [1, 8, 25, 116, 119], this coverage can be traded off to 
focus on some specific regions, in favor of increased func- 
tional specificity and thus better differentiation across net- 
works [28, 46, 114]. In addition, while ROIs are most often 
defined as a hard selection of voxels, it is also possible to 
use a soft definition, attributing weights as with proba- 
bilistic atlases, or spatial maps of functional networks ex- 
tracted from techniques such as independent component 
analysis (ICA) [57, 99]. 

Regions from atlases. Atlases can be used to define full- 
brain par cellat ions. Popular choices are the Automatic 
Anatomic Labeling (AAL) atlas [111], which benefits from 



1 In neuroimaging, the term network is sometimes used to denote 
a graph of brain function. To disambiguate the notion of segregated 
spatial mode [31] from that of connectivity graphs, we will purposely 
restrict its usage in this paper. 
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Note that naive use of filtering can induce spurious correlations 



an SPM toolbox, or the ubiquitous Talaraich-Tournoux at- 
las [107]. However, these atlases suffer from major short- 
comings; namely i) they were defined on a single subject 
and thus do not reflect inter-subject variability, and ii) 
they focus on labeling large anatomical structures and do 
not match functional layout -for instance only two re- 
gions describe the medial part of the frontal lobe in the 
AAL atlas. Multi-subject probabilistic altases such as 
the Harvard- Oxford atlas distributed with FSL [98] or the 
sulci-based structural atlas used in [116] mitigate the first 
problem, and the high number of regions defined using 
sulci also somewhat circumvent the second problem (see 

fig- 1)- 

Defining regions from the literature. Regions can be de- 
fined from previous studies, informally or with system- 
atic meta-analysis. This strategy is used to define the 
main resting- state networks, such as the default mode net- 
work, but may also be useful to study connectivity in task- 
specific networks [14, 28, 47, 86]. The common practice is 
to place balls of a given radius, 5 or 10 mm, centered at the 
coordinates of interest. Given that functional networks are 
tightly interleaved in some parts of the cortex, such as the 
parietal lobe, care must be taken not to define too many 
regions that would overlap and lead to mixing of the signal. 

FMRI-based function definition. Defining regions directly 
from the fMRI signal brings many benefits. First, it can 
capture subject-specific functional information. Second, 
it adapts to the signal at hand and its limitations, such 
as image distortions or vascular and movement artifacts 
that are isolated in ICA-like approaches. Lastly, incorpo- 
rating functional information into regional definition will 
result in more homogenous regions that better represent 
connectivity present at the voxel level than anatomically- 
defined atlases such as AAL or Harvard- Oxford [25]. The 
simplest approach to define task-specific regions is to use 
activation maps derived from standard GLM-based anal- 
ysis in a task-driven study (see for instance [81]). Re- 
gions are extracted by thresholding the maps, or using 
balls around the activation peaks. For resting- state stud- 
ies, unsupervised multivariate analysis techniques are nec- 
essary. Clustering approaches extract full-brain parcella- 
tions [9, 25, 109, 121], and have been shown to segment 
well-known functional structures from rest data. Alterna- 
tively, decomposition methods, such as ICA [6], can unmix 
linear combinations of multiple effects and separate out 
partially-overlapping spatial maps that capture functional 
networks or confounding effects, as for instance with the 
presence of vascular structure in functional networks. At 
high model order, ICA maps define a functional panel- 
lation [57]. Extracting regions from these maps requires 
additional effort as they can display fragmented spatial 
features and structured background noise, but incorporat- 
ing sparsity and spatial constraints in the decomposition 
techniques leads to contrasted maps that outline many dif- 
ferent structures [117] (see fig. 1). 
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Figure 1: Different full-brain parcellations: the AAL atlas [111], the 
Harvard- Oxford atlas, the sulci atlas used in [116], regions extracted 
by Ncuts [25], the resting-state networks extracted in [97] by ICA, 
and in [115] by sparse dictionary learning. 



Optimal number of regions. Defining an optimal number 
of regions to use for whole-brain connectivity analysis 
bears careful consideration. On one hand we desire a suf- 
ficiently large number of regions to guarantee that they 
are functionally homogeneous regions and adequately rep- 
resent the connectivity information present in the data. 
On the other hand too many regions will render statistical 
inference challenging, result in an explosion in computa- 
tional complexity, and interfere with the interpretability of 
observed connections. For functional parcellation, cross- 
validation methods can be employed to estimate an opti- 
mal number of regions based on homogeneity, the ability 
to reproduce connectivity information present at the voxel 
scale, and the ability to obtain the same parcellations from 
independent data [15, 25]. In general these metrics do not 
result in an obvious peak at a "best" number of regions, 
but instead offer a range over which the number of regions 
can be chosen based on the needs of the analysis at hand. 
Finally, it is important to keep in mind that there is no 
universally better parcellation and associated number of 
regions. From a practical standpoint, these choices will 
depend on the task at hand, and more fundamentally, a 
good description of brain function should cover multiple 
scales. Given that it is not clear that an optimal parcel- 
lation can be identified from the sample size of a typical 
study, randomized parcellation, as used in structural con- 
nectomes [124] or activation mapping [118], may also be 
considered. 

2.3. Estimating connections 

The concept of functional connectivity has been called 
elusive [51]: it has many mathematical instantiations al- 





Figure 2: Correlation matrices of rest time-series extracted from the 
39 main regions of the Varoquaux 2011 [115] parcellation with differ- 
ent choices of confound regressors - Left: regressing out CompCor 
signals, as well as white matter and CSF average signals and move- 
ment parameters. The insert shows the connections restricted to a 
few major nodes. - Upper right: regressing out only movement pa- 
rameters. - Lower right: regressing out movement parameters and 
global signal mean. No frequency filtering was applied here. When 
no confounding brain signals are regressed, all regions are heavily cor- 
related. Regressing out common signal, in the form of well- identified 
confounds or a global mean, teases out the structure. 



though in essence they all strive to extract simple statistics 
from functional imaging in order to characterize synchrony 
and communication between large ensembles of neurons. 
Here we choose to focus on second order statistics that 
can be related to Gaussian models, the simplest of which 
being the correlation matrix of the signals of the different 
ROIs. 

Signal extraction. Given a set of graph nodes, the next 
step is to extract a representative time series for each node. 
To study intrinsic activity, e.g. with rest data, signal ex- 
traction can be achieved by either averaging the fMRI time 
series across the voxels in a region, or by taking the first 
eigenvariate from a principle components analysis of the 
time series [40] . Comparisons of these methods has shown 
that the eigenvariate method is more sensitive to function 
inhomogeneity [25] and exhibits worse test-retest reliabil- 
ity than averaging time series [128]. In addition, improved 
specificity to BOLD signal can be enforced by using only 
signal in voxels near gray-matter tissues. For this pur- 
pose, we suggest summarizing the signal in an ROI by 
a mean of the different voxels weighted by the subject- 
specific gray matter probabilistic segmentation, as output 
by e.g. SPM's segmentation tool [4] or FSL's FAST pro- 
gram [126]. 

Studying connectivity from evoked activity with task- 
driven studies requires disambiguating task-specific con- 
nectivity effects from intrinsic connectivity mediated by 
shared neuromodulatory/task inputs, anatomical path- 
ways, etc. In this regard, it can be beneficial to run a 
GLM-based first-level analysis, enforcing specificity of the 
measure extracted to the task. With slow event-related 



Figure 3: Different inverse-covariance matrices estimates corre- 
sponding to fig. 2 - Left: group-sparse estimate using the &21 es- 
timator [116]. The insert shows the connections restricted to a few 
major nodes. - Upper right: non-sparse estimate: inverse of the 
sample correlation matrix. - Lower right: sparse estimate using 
the Graph Lasso [34]. 



designs, task-specific functional connectivity can be cap- 
tured in trial-to-trial fluctuations in the BOLD response, 
estimated using a GLM analysis with one regressor per 
trial [47, 75, 86]. This approach, known as beta-series re- 
gression, has been adapted for rapid event-related designs, 
using multiple GLMs to optimize deconvolution of each 
trial [76]. 

Correlation and partial correlations. Given ROIs defining 
the nodes of the functional-connectome graph, one needs 
to estimate the corresponding edges connecting them. 
Functional connectivity between the ROIs can be mea- 
sured by computing the correlation matrix of the extracted 
signals. An important and often neglected point is that the 
sample correlation matrix, i.e. the correlation matrix ob- 
tained by plugging the observed signal in the correlation 
matrix formula, is not the population correlation matrix, 
i.e. the correlation matrix of the data-generating process. 
If the number of measurements was infinite, the two would 
coincide, however if this number is not large compared to 
the number of connections (that scales as the square of the 
number of ROIs), the sample correlation matrix is a poor 
estimate of the underlying population correlation matrix. 
In other words, the sample correlation matrix captures a 
lot of sampling noise, intrinsic randomness that arises in 
the estimation of correlations from short time series. Con- 
clusions drawn from the sample correlation matrix can eas- 
ily reflect this estimation error. Varoquaux et al. [116] and 
Smith et al. [100] have shown respectively on rest fMRI and 
on realistic simulations that a good choice of correlation 
matrix estimator could recover the connectivity structure, 
where the sample correlation matrix would fail. In general, 
the choice of a better estimate depends on the settings and 
the end goals [114, 117], however the Ledoit-Wolf shrink- 
age estimate [62] is a simple, computationally-efficient, and 




Figure 4: Inverse-covariance matrices for different choice of con- 
found regressors - Left: regressing out only movement parameters 
- Right: removal of the global mean, instead of the white matter, 
CSF, and CompCor time courses. 



parameter-less alternative that performs uniformly better 
than the sample correlation matrix [116, 117] and should 
always be preferred. 

For the problem of recovering the functional- 
connectivity structure, i.e. finding which region is con- 
nected to which, sparse inverse covariance estimators have 
been found to be efficient [89, 100, 116]. The intuition for 
relying on inverse covariance rather than correlation stems 
from that fact that standard correlation (marginal correla- 
tion) between two variables a and b also capture the effects 
of other variables: strong correlation of a and b with a third 
variable c will induce a correlation between a and b. On 
the opposite, the inverse covariance 3 matrix (also called 
precision matrix) captures partial correlations, removing 
the effect of other variables [71]. In the small sample limit, 
this removal is challenging from the statistical standpoint. 
This is why an assumption of sparsity, i.e. that only few 
variables need to be considered at a time, is important 
to estimate a good inverse covariance. Various estimation 
strategies exist for sparse inverse covariance, and have an 
impact on the resulting networks [116, 117]. The Graph 
Lasso (£i -penalized maximum-likelihood estimator) [34] is 
in general a good approach for structure recovery. In group 
studies, the £21 estimator [50, 116] is useful to impose a 
common sparsity structure across different subjects and 
achieve better recovery of this common structure. Simply 
put, these approaches are necessary because estimation 
noise creates a background structure (see fig. 3); however, 
unlike in a univariate situation, the parameters are not in- 
dependent, and the spurious background connections de- 
grade the estimation of the actual connections. The sparse 
estimators make a compromise between imposing simpler 
models, i.e. with less connections, and providing a good 
fit to the data. This compromise is set via a regulariza- 
tion parameter which controls the sparsity of the estimate. 
A good procedure to choose this parameter is via cross- 



3 Covariance and correlation matrices differ simply by the fact 
that a covariance matrix captures the amplitude of a signal, via its 
variance, while a correlation matrix is computed on standardized 
(zero mean, unit variance) signals. 



validation [116]. 

Network structure extracted. The correlation matrices and 
inverse-covariance matrices that we extract contain a lot of 
information on the functional structure of the brain. First, 
the correlation matrix (fig. 2) shows blocks of synchronized 
regions that can be interpreted as large-scale functional 
networks, such as the default mode network. Note that the 
split in networks is not straightforward. Different ordering 
of the nodes will reveal different networks. Indeed, because 
of the presence of hubs and interleaved networks, the pic- 
ture in terms of segregated networks is not sufficient to ex- 
plain full-brain connectivity [117]. Connectivity matrices, 
correlation matrices and inverse-covariance matrices, can 
be represented as graphs: nodes connected by weighted 
edges (inserts on fig. 2 and fig. 3). The inverse-covariance 
matrix, which captures partial correlations, appears then 
as extracting a backbone or core of the graph. While such 
structure has been used as a way to summarize anatomi- 
cal brain connectivity graphs [49], here it has a clear-cut 
meaning with regard to the BOLD signal: it gives the con- 
ditional independence structure between regions [117]. In 
other words, regions a and b are not connected if the sig- 
nal that they have in common can be explained by a third 
region c. In this light, the choice of nuisance regressors 
to remove confounding common signal is less critical with 
partial correlations than with correlations. Indeed, while 
with correlation matrices regressing out the global mean 
has a drastic effect (fig. 2 upper right and lower right), on 
inverse covariance it only changes the resulting matrices 
very slightly (fig. 4). 

There have been debates on whether to regress out cer- 
tain signals, such as the global mean, as it induces negative 
correlations [19, 32, 77], and these may seem surprising: 
one network appears as having opposite fluctuations to 
another. However, correlation between two signals only 
takes its meaning with the definition of a baseline. A sim- 
ple picture to explain anti-correlations between two regions 
is the presence of a third region, mediating the interac- 
tions. Using this third region as a baseline would amount 
to estimate partial correlations in the whole system. Us- 
ing inverse-covariance matrices or partial correlations to 
understand brain connectivity makes the interpretation 
in terms of interactions between brain regions easier and 
more robust to the choice of confounds. 



3. Comparing connectivity 

We now turn to the problem of comparing functional 
connectivity across subjects or across conditions. 

3.1. Detecting changing connections 

First, we focus on detecting where the connectivity ma- 
trices estimated in the previous section differ. 



Mass-univariate approaches. The most natural approach 
is to apply a linear model to each coefficient of the con- 
nectivity matrices [47, 64]. This approach is similar to the 
second-level analysis used in mass-univariate brain map- 
ping, and gives rise to many of the well-known techniques 
used in such a context, such as the definition of a second- 
level design, with possibly the inclusion of confounding ef- 
fects, and statistical tests (T tests or F tests) on contrast 
vectors. Importantly, in order to work with Gaussian- 
distributed variables, it is necessary to apply a Fisher Z 
transform 4 to the correlations. Note that in these set- 
tings, the Ledoit- Wolf estimator [62] is often a good choice 
to estimate the correlation matrix, as it is parameter-free 
and gives good estimation performance without imposing 
any restrictions on the data. For hypothesis testing, cor- 
recting for multiple comparison can severely limit statis- 
tical power, as the number of tests performed scales as 
the square of the number of regions used. Controlling for 
the false discovery rate (FDR) mitigates this problem. Al- 
ternatively, as the assumptions underlying the Benjamini- 
Hochberg procedure [10] for the FDR can easily be broken, 
non-parametric permutations-based tests give reliable ap- 
proaches. In particular, the max-T procedure [42, 79] is 
interesting to avoid the drastic Bonferroni correction when 
controlling for multiple comparison in family-wise error 
rate. 

Accounting for distributed variability. A specific challenge 
of connectivity analysis is that the connectivity strength 
between different regions tends to covary. For instance, 
with resting-state data, functional networks comprising 
many nodes can appear as more or less connected across 
subjects (see for instance fig. 5, showing variability in a 
control population at rest). In other words, non-specific 
variability is distributed across the connectivity graph, 
and it is structured by the graph itself. This obser- 
vation brings the natural question of whether second- 
level analysis should be performed on correlation matri- 
ces, inverse- covariance matrices, or another parametriza- 
tion that would disentangle effects and give unstructured 
(white) residuals. While inverse-covariance matrices show 
less distributed fluctuations than correlation matrices, 
they capture a lot of background noise, as partial corre- 
lations are intrinsically harder to estimate. Preliminary 
work [114] suggests performing statistical tests on residu- 
als of a parametrization intermediate between correlation 
matrix and inverse covariance matrix, as it can decouple 
effects and noise. 

Taking a different stance on distributed variability, the 
"network-based statistics" approach [122] draws from the 
hypothesis that if, in a second- level analysis, an effect is 
detected on a connection that lies in a network of strongly 
connected nodes, a large sub-network is likely to carry an 




Figure 5: Inter subject variability. Note that this is variability occur- 
ring in a healthy population at rest, in other words it is non specific 
variability - a: single-subject correlation matrices for different sub- 
jects - b: Corresponding Z-score (effect / standard deviation) of the 
difference between a subject and the remaining others - c: single- 
subject inverse-covariance matrices - d: Corresponding Z-score for 
the inverse-covariance matrices - e: Corresponding Z-score for the 
subject residuals, as denned in [114]. 



effect. Thus, they adapt cluster- level inference to connec- 
tivity analysis, in order to mitigate the curse of multiple 
comparisons. 

3.2. Comparing network summary statistics 

Both the multiple comparison issue and the network- 
level distributed variability are a plague to edge-level com- 
parison of connectomes. A possible strategy to circumvent 
these difficulties is to perform comparisons and statistical 
testing at the level of the network, rather than the indi- 
vidual connection. 

Network integration. Marrelec et al. [69] introduce the 
use of entropy and mutual information as a measure of 
network- level functional integration 5 . Gaussian entropy 
can be seen as a simple metric to generalize correlation or 
variance to multiple nodes (see [3] §2.5.2 and §7.5). In- 
deed, let us consider 3 nodes: a, b and c. Their correlation 
structure is captured by three correlation coefficients: p a 5, 
pb c and p ac . Summarizing these by their mean, as might 
seem natural, discards the relationship between the sig- 
nals, while using the integration metric, defined as the 
Gaussian entropy, tells us how much two signals can be 



4 See http://en.wikipedia.org/wiki/Fisher_transformation or 
[3] section 4.2.3 for mathematical arguments. 



5 See [116] for simplified formulas for network integration and mu- 
tual information. 
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0.554 




Integration: 
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Figure 6: Two different correlation matrices with the same average 
correlation, but with very different integration values. Indeed, the 
matrix on the left was chosen to represent three signals a, b and c as 
different from each other as possible, given p a ^ + p^ c + p ac = 1.35; it 
thus has a small integration value. On the opposite, for the matrix 
on the right, signal b can almost be fully recovered by combining 
signals a and c; the matrix thus has a large integration value. 



combined to form the third (see fig. 6). Cross-entropy -or 
mutual information- [69] measures the amount of cross- 
talk between two systems in a similar way as Gaussian 
entropy is used to measure the integration of a brain sys- 
tem. The functional-connectivity structure, or its repre- 
sentation in the form of a correlation matrix, can thus be 
characterized via the integration and cross-talk of some 
of its sub- systems. This approach gives a simplified rep- 
resentation with a small number of metrics that can be 
compared across subjects. 

Graph-topological metrics. Functional connectivity graphs 
have been found to display specific topological prop- 
erties 6 that are characteristic of small-world networks 
[1, 17, 91, 103]. These networks display excellent transport 
properties: although they have a relatively small number 
of connections, any two regions of the brain are well con- 
nected. Another interesting consequence of their specific 
topology is the resilience it gives the system to attacks such 
as resulting from brain lesions [1]. This overall structure 
of functional-connectivity graphs can be summarized by 
a few metrics, such as the average path length between 
any two nodes, the local clustering coefficients, or the 
node degree centrality [87]. Given that pathologies with- 
out a localized focus, such as schizophrenia, are thought 
to have a global impact on brain connectivity [5, 65], the 
graph-topological metrics are promising markers to per- 
form inter-subject comparison. Such an approach is ap- 
pealing as it is not subject to multiple comparison issues. 
However, it has been criticized as giving a fairly unspe- 
cific characterization of the brain and being fragile to noise 
[54]. Another caveat is that these properties are not spe- 
cific to brain function: correlation matrices display small- 
world properties such as local clustering by construction. 
Indeed, if two nodes are strongly correlated to a third, 
they are highly likely to be correlated to each other [123]. 
This observation highlights the need for well defined null- 
hypothesis [88, 123], but also for controlled recovery of 
brain functional connectivity going beyond empirical cor- 
relation matrices, as discussed in the previous section. 



3.3. Predictive Modeling 

Predictive modeling is concerned with learning (or fit- 
ting) a model that is capable of predicting information 
from unseen data [80]. In the context of connectomes, pre- 
dictive modeling can extract connectivity-based biomark- 
ers of disease diagnosis, prognosis, or other phenotypic 
outcomes [24, 27]. The accuracy of a predictive model 
provides a measure of the amount of information present 
in the connectome about the phenotypic measure being 
evaluated [58, 59]. When combined with reproducibility, 
prediction accuracy provides a metric for evaluating ex- 
perimental trade-offs for data acquisition, preprocessing, 
and analysis [60, 106]. Multivariate predictive models are 
attractive in connectomics because they are sensitive to 
dependencies between features and avoid the need to cor- 
rect for multiple comparisons since the significance of an 
entire pattern is evaluated using a single statistical test. 
Additionally, modern predictive modeling techniques draw 
from the statistical learning literature, which specifically 
addresses high dimensional datasets with few observations. 
Predictive modeling has been successfully applied to iden- 
tify connectome-based biomarkers of Alzheimer's disease 
[104], depression [24, 125], schizophrenia [18, 94], autism 
[2], ADHD [127], aging [27], as well as to classify mental 
operations [85, 95]. The growing interest for applying pre- 
dictive modeling to connectivity analysis was highlighted 
by the ADHD200 Global Competition, in which the object 
was to identify a connectivity-based biomarker of ADHD 
[108]. Recent work has illustrated the utility of predictive 
modeling for deriving connectivity models at the individ- 
ual level [22]. 

Technically, predictive modeling is a supervised ma- 
chine learning problem where a target to be predicted 
-e.g. age, disease state, cognitive state- is available for 
each observation of the data. In the context of comparing 
connectomes, features used in the predictive model corre- 
spond to bivariate measures of connectivity [27, 85, 95], 
or any of the previously discussed graph summary metrics 
[18, 29]. The quality of a predictive model is determined by 
its prediction accuracy (or generalization ability) which is 
measured using one or more iterations of cross-validation. 
Cross-validation iteratively subdivides available data into 
a subset used for training the classifier and a dataset for 
evaluating classifier performance 7 [80] . The significance of 
achieved prediction accuracy can be assessed using permu- 
tation tests [44] . Predictive modeling approaches typically 
require the specification of several parameters, which may 
be chosen based on domain specific knowledge or require- 
ments [21], determined using an analytical approach [20], 
or optimized using a second-level cross validation proce- 
dure [33]. 



6 In the neuroscience world, these descriptions are grouped under 
the terms of "graph-theoretical approaches", however graph theory 
is an entire division of mathematics and computer science that is 
concerned with much more than topology of random graphs. 



7 Several strategies exist for performing cross-validation and the 
commonly used approach of using only a single observation for testing 
(leave-one-out cross-validation) results in highly variable estimates 
of prediction accuracy [33]. Alternative approaches such as (5 or 10)- 
fold cross-validation, or 0.632+ bootstrap should be preferred [33]. 



Although predictive modeling techniques are well 
suited for measuring whether information exists in the con- 
nectome about a phenotypic variable, they do not directly 
identify the connections that are most relevant to the pre- 
diction. This limitation can be somewhat mitigated by 
relying on previously described sparse inverse covariance 
estimation techniques to minimize the number of connec- 
tions. Additionally, feature selection [48] can be performed 
by filtering out features based on their statistical relation- 
ship with the variable of interest [24, 94]. Importantly, it 
must be performed within cross-validation to avoid bias- 
ing estimates of prediction accuracy. The interpretation of 
connections used in predictive models and their relation- 
ship with a phenotypic outcome is difficult and requires 
insight into the mechanisms underlying the modeling ap- 
proach. For linear models, the weights of the model are 
similar to the weights of an ordinary linear regression. If 
features are appropriately standardized prior to training, 
the magnitude of the weights can be interpreted to re- 
flect the relative importance of the feature to the model, 
for instance the corresponding edge in the connectome. 
However, interpreting how the connections differ between 
classes or relate to a phenotypic variable can be more com- 
plicated given the multivariate nature of their involvement. 
Indeed, the inclusion or exclusion of a connection in the 
model can induce a change of the sign of another model 
weight [24] . It is perhaps most reasonable to adapt a con- 
servative interpretation in which predictive modeling is 
used to identify candidate connections that are later tested 
in follow-up experiments better suited to elucidating their 
relationship to the variable of interest. 

To conclude on predictive modeling with a practical 
note for connectome comparison, we would like to stress 
that while machine learning algorithms are powerful tools, 
they work best if they are provided with discriminant- 
ing and noiseless features. In other words, as with all 
other connectome comparison methods, optimizing first- 
level analysis -subject-level connectome extraction- is 
paramount. 

4. Beyond correlation, effective connectivity? 

All the approaches that we have presented in this re- 
view are based on second-order statistics of the signal, 
in other words correlation analysis. Traditionally, these 
are defined as functional connectivity, defined as "tempo- 
ral correlations between remote neurophysiological events" 
[35], and opposed to effective connectivity, i.e. "the in- 
fluence one neural system exerts over another" [35]. To 
conclude this review, we would like to bridge the gap be- 
tween these concepts, which in our eyes should be seen as 
a continuum rather than an opposition (this opinion is also 
expressed in [73]). 

A first step to move from purely descriptive statistics 
to interaction models with functional connectivity analysis 
is to consider a correlation matrix as a Gaussian graphi- 
cal model, i.e. a well-defined probabilistic model that de- 



scribes observed correlations in terms of an independence 
structure and conditional relations [61, 117]. In such set- 
tings, the inverse covariance graph or the partial correla- 
tions are a measure of influence from one node to another, 
albeit undirected. Inferring directionality in a Gaussian 
model is impossible. Linear structural equation models 
(SEMs) [74] rely on a similar model that consists in speci- 
fying a candidate directed graphical structure. This struc- 
ture constraints the covariance matrix of the signals and 
can thus be tested on observed data. In fact some forms of 
SEMs are known as "covariance structure models" . There 
is thus a strong formal link between correlation analysis in 
the framework of graphical models and SEMs: the former 
is undirected but fully exploratory, as it does not require 
the specification of candidate structure, while the latter is 
directed but confirmatory. This link has been exploited to 
specify candidate structures for SEMs using partial cor- 
relations [70]. More complex models, such as dynamical 
causal models (DCMs) [38] or Granger causality [43] re- 
quire additional hypotheses such as non-linear couplings 
or time lags. 

Most importantly, more complex models can only be 
used to model interactions between a small number of 
nodes. This is not only due to a computational difficulty, 
but also to fundamental roadblocks in statistics: the com- 
plexity of the model must match the richness of the data. 
While injecting prior information can help model estima- 
tion, the more informative this prior is, the more fragile 
the inference becomes. The ongoing debate on the impact 
of hemodynamic lag on Granger-causality inference [96] 
is an example of such fragility. Note that although most 
of the theory underpinning correlation analysis (Gaussian 
graphical models) is based on a Gaussian assumption, the 
core results are robust to violations of this assumption [84] . 

It is tempting to favor more neurobiologically-inspired 
models that give descriptions close to our knowledge of the 
brain basic mechanisms, however, as George Box famously 
said, "all models are wrong; some models are useful" . De- 
pending on the question and the data at hand, a trade-off 
should be chosen between complex models based on a bio- 
physical description, and simple phenomenological models 
such as correlation matrices. In particular to model inter- 
actions between a large number of regions, as in full-brain 
analysis, and learn a large connectome, simple models are 
to be preferred. For more hypothesis-driven studies, such 
as the analysis of the mechanisms underlying a specific 
task, more complex models can be preferred, if rich data 
is available. Automatic choice of model is a difficult prob- 
lem, however, cross-validation (as used in [25, 105, 116]) 
is a useful tool. The central principle of cross-validation is 
to test a model on different data than the data used to fit 
the model. Models too complex for the data available will 
fit noise in the data, and thus generalize poorly. The main 
benefit of cross-validation is that it is a non-parametric 
method which does not rely strongly on modeling assump- 



tions . 



5. Conclusion 

Horwitz el al. [52] claimed almost 20 years ago that 
"the crucial concept needed for network analysis is covari- 
ance" . In our eyes, this still holds today. Estimation func- 
tional connectomes relies largely on fitting covariance mod- 
els. Their comparison requires understanding how these 
covariances vary and finding metrics to capture this vari- 
ability. The additional secret ingredient may be using con- 
founds regressors in all statistical steps. A good choice of 
a small number of relevant regions facilitates connectome 
comparison. However, such a choice cannot yet be fully 
factored out via methods and must rely on neuroscientiflc 
expertise. 

Methodological challenges to functional-connectome- 
based group studies arise from the dimensionality and 
the variability of the connectome. With the current 
tools, inter-subject comparison of connectomes compris- 
ing many nodes is limited by the difficulty of estimating 
high-dimensional covariance matrices and the loss of sta- 
tistical power due to multiple comparisons. Better algo- 
rithms integrating powerful a priori information are re- 
quired to push the limits of covariance estimation. Better 
characterization of inter-subject variability of connectomes 
[56] will help choosing parameterizations and invariants to 
avoid testing each edge for a difference, as this strategy 
inevitably leads to a needle in a haystack problem. 

Reviewing methodological options to learn and com- 
pare connectomes highlights that there is currently no 
unique solution, but a spectrum of related methods and 
analytical strategies. More empirical results are required 
to guide the choices. However this diversity is probably 
unavoidable: a diffuse disease like schizophrenia will not 
lead to the same connectome modifications as a focal le- 
sion. In statistical learning, "no free lunch" theorems [120] 
tell us that no strategy can perform uniformly better in all 
situations. In practice, the key to a successful analysis is 
to understand well the assumptions and interpretation of 
each option, in order to match the method to the question. 
Similarly, the idealized notion of an unique functional con- 
nectome to describe connections in brain function is prob- 
ably an utopia, and various connectomes should be con- 
sidered in different settings, such as the study of varying 



8 This is to be contrasted to Bayesian model comparison, which 
will give well-controlled results only if the true generative model is in 
the list of models compared. [36] argues that, based on the Neyman- 
Pearson lemma, cross-validation is less powerful than likelihood ratio 
tests using the full dataset. However, it is important to keep in mind 
that these approaches only test for self-consistence, as the Neyman- 
Pearson lemma is established under the hypothesis that the model 
used to define the test is indeed the data-generating process [78], 
while in practice it is often the case that this model gives poor fits to 
the data [66]. Applying test procedures on different data than that 
used to fit the model, as in cross-validation, is much more resilient 
to modeling errors. 



phenotypic conditions, or that of on-going activity versus 
activity related to specific tasks. 
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